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Abstract 

A recently developed method, introduced in Phys. Rev. Lett. 94 (2005) 180403, Phys. 
Rev. B 72 (2005) 064302, Phys. Lett. A 344 (2005) 84, systematically improved the 

convergence of generic path integrals for transition amplitudes. This was achieved 

(p) 

by analytically constructing a hierarchy of iV-fold discretized effective actions S*r 
labeled by a whole number p and starting at p = 1 from the naively discretized action 
in the mid-point prescription. The derivation guaranteed that the level p effective 
actions lead to discretized transition amplitudes differing from the continuum limit 
by a term of order 1/N P . Here we extend the applicability of the above method to 
the calculation of energy expectation values. This is done by constructing analytical 
expressions for energy estimators of a general theory for each level p. As a result 
of this energy expectation values converge to the continuum as 1/N P . Finally, we 
perform a series of Monte Carlo simulations of several models, show explicitly the 
derived increase in convergence, and the ensuing speedup in numerical calculation 
of energy expectation values of many orders of magnitude. 
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1 Introduction 



The path integral formalism first introduced by Feynman [4,5] represents a 
rich and flexible general mathematical setting for dealing with quantum and 
statistical theories. The most obvious success of path integrals has been the 
ease with which they have allowed us to extend the quantization procedure 
to ever more complicated systems. On a more practical note, the formalism 
has been extremely useful for handling symmetries, deriving non-perturbative 
results and connections between different theories [6,7], and as a catalysts for 
the exchange of key ideas between different areas of physics, most notably high 
energy and condensed matter physics [8,9]. Today, analytical and numerical 
approaches to path integrals [10,11,12,13] play important roles not only in 
physics but also in chemistry, materials science, mathematics and modern 
finance. 

Further development of the path integral method is constrained by the small 
number of solvable models, as well as by our rather limited knowledge of their 
precise mathematical properties. In an attempt to fill this void a recent series 
of papers [1,2] has investigated the dynamical implications of the property of 
stochastic self-similarity of path integrals by studying the relation between 
path integral discretizations of different coarseness. This has resulted in a 
systematic analytical construction of a hierarchy of iV-fold discretized effec- 
tive actions S$ labeled by a whole number p and built up from the naively 
discretized action in the mid-point prescription (corresponding to p — 1). It 
was shown that the level p effective actions lead to discretized transition am- 
plitudes differing from the continuum limit by a term of order 1/N P . These 
analytical results in fact represent an explicit derivation of the generalization 
of Euler's summation formula to path integrals [3]. 

From a numerical stand point the new method has lead to a many order of 
magnitude speedup of path integral simulations by making it possible to get 
precise results using small values of N. The substantial numerical speedup is 
the direct result of the new analytical input that makes it possible to trade 
high N (high cost in computing time) for low iV and high p. For example, for 
p = 9 the overall speedup in convergence of the algorithm is typically eight 
orders of magnitude [2] over the defining algorithm, or five orders of magnitude 
when compared to previous state-of-the-art algorithms [14,15,16] based on the 
generalized Trotter formula [17], or other short-time approximation schemes 
[18]. 

In this paper we extend the applicability of the above method from transition 
amplitudes to the calculation of expectation values. In particular, we focus 
on energy expectation values. It is well known that the efficient calculation of 
energy expectation values necessitates the use of the optimal energy estima- 
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tors [19,20,21,22]. An excellent review of different energy estimators is given in 
[23] . In order to increase convergence of energy expectation values to the con- 
tinuum limit it is necessary to understand how the estimators need to change 
as we move through the hierarchy of effective actions. In fact, if estimator 
and effective action are not synchronized then the increase in convergence for 
transition amplitudes does not translate into a corresponding increase in con- 
vergence for energy expectation values. The central result of this paper is the 
explicit analytical construction of optimal energy estimators for each hierarchy 
level p. We show that energy expectation values calculated using p level effec- 
tive actions and the associated energy estimators derived here converge to the 
correct continuum value as 1/N P , i.e. lead to the same increase in convergence 
(and algorithm speedup) previously seen for transition amplitudes. We have 
carried out a series of Monte Carlo simulations of several models (anharmonic 
oscillator, Poschl- Teller potential, Morse potential) and have shown explicitly 
that the increase in convergence (and the ensuing speedup in numerical cal- 
culation of energy expectation values) agrees with the analytical derivations. 
The figures presented in the paper illustrate the results of these simulations 
for the case of an anharmonic oscillator with quartic coupling. 

Section 2 starts with a brief overview of notation and an introduction to 
different energy estimators. It then develops the procedure for determining 
the estimator associated to each effective action at hierarchy level p. Section 3 
presents the results of our Monte Carlo simulations. In the Appendix we give 
the explicit expressions for the energy estimators for p < 6. The results for 
p < 9 and the codes used can be found on our web site [24] . 



2 Energy estimators 



In the functional formalism the transition amplitude (in Euclidean time) 
A(a, b; T) = {b\e~ TH \a) is given in terms of a path integral which is simply the 
N — > oo limit of the (N — 1)- fold integral expression 

1 ~ 

A N (a, b; T) = (^^) * / d qi ■ ■ ■ dq N ^ e~ SN . (1) 



The Euclidean time interval [0, T] has been subdivided into N equal time steps 
of length €n = T/N, with qo = a and qN = b. The integrand is given in terms 
of the naively discretized action Sn- For actions of the form 

T 

S = jdt (^q 2 + V{q)) , (2) 
o 
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the naively discretized action (in the mid-point ordering prescription) equals 

N-l / p \ 

S n = E + tNV(q n )j , (3) 

where S n = q n+ \ — q n , and q n = |(? n +i + ?n)- Note that we use units in which 
h and particle mass have been set to unity. 

As can be seen, the very definition of path integrals makes necessary the tran- 
sition from continuum to the discretized theory. This discretization, however, 
is far from unique. We are free to introduce additional terms that explic- 
itly vanish in the continuum limit. Although such additional terms do not 
change the continuum physics, they do affect the speed of convergence to that 
continuum limit. A recent series of papers [1,2] has studied the relation be- 
tween discretizations of different coarseness and has analytically constructed 
a hierarchy of effective actions S$. The starting member of the hierarchy 
corresponds to p — 1 and is given by Eq. (3), the naively discretized action 
in the mid-point prescription. The p-level effective action leads to discretized 
amplitudes that converge to their continuum values as 

A® (a, b; T) = A(a, b;T) + {{T/Nf) . (4) 



Before extending the above outlined procedure to improve the convergence of 
energy expectation values, we briefly review the standard procedure for their 
calculation. The canonical partition function of statistical mechanics 

Z(/3) = Tre^ = J dq A(q,q; (3) , (5) 



is given in terms of diagonal transition amplitudes where the inverse tempera- 
ture f3 plays the role of propagation time T. From the above it directly follows 
that the partition function may be written in terms of a path integral - the 
continuum limit of the TV-fold integral 

1 ~ 

Z N (P) = (7^—) 2 fd qi --- dqN^dqN e~ SN , (6) 

where integration is over all periodic trajectories q(t) = q(t + 13). In the dis- 
cretized theory this simply implies that q = q N . 

Path integral Monte Carlo simulations use two different types of energy esti- 
mators - "kinetic" and "virial" . Both types of estimators are well known and 
have been extensively studied in [19,20,21,22,23]. We start with the "kinetic" 
estimator which follows from straightforward differentiation of the partition 
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function. The (thermal) expectation value of the energy U{(3) is simply the 
continuum limit of 



U N (J3) = ~\nZ N (j3) . 



(7) 



Note that (3 = New, so that d/d/3 = (l/N)d/deN- Using this we finally find 
that 



Z]y ^\2lT€N 



dqi... q N 



1 dS N 
+ 



2€n de 



N 



Putting in the naively discretized action we recover the so-called "kinetic" 
estimator of the energy 



N 



N-l 



N-l 



-^kinetic ~ IT AT ^ + AT ^ ^^ n ^ ' 



N 



(9) 



n=0 



The last term on the right hand side is clearly the potential energy estimator. 
The second term is what one would naively take to be the estimator for kinetic 
energy. This term, however, diverges in the continuum limit. In fact, the role of 
the first term in the above expression (the one coming from the path integral 
measure) is precisely to cancel this divergence. The above estimator, taken as a 
whole, is well behaved in the continuum limit. However, the fact that it is given 
as a difference of two diverging terms implies that its variance diverges with 
N. It is advantageous, therefore, to find another estimator having the same 
mean value but with smaller variance. This second type of energy estimator 
is called "virial" as it is based on the virial theorem 



(xV'(x)) 



(10) 



To derive this in the path integral formalism we rescale the coordinates in 
Eq. (6) according to q n — > /iq n . By setting /i 2 = all the f3 dependence is 
removed from the path integral measure, and we get 



N-l 



z n(P) = Q^j J dq ± ... dq N exp j- ^ 



5, 



f + e N V(e 2 N q n ) 



;n) 



Differentiating with respect to we find the "virial" estimator for the energy 
to be 

i N-l i N-l 

^virial = OAF E ^V'(qn) + - £ V(q n ) . (12) 

ZiV n=0 iV n=0 
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As expected, the new estimator is a sum of terms that are finite in the con- 
tinuum limit and so leads to smaller variance. The remaining work presented 
in this paper uses the "virial" estimator given above and its p level general- 
izations. From now on we drop the designation "virial" . 

As we have seen, the estimator follows directly from the naive discretized 
action. If instead of Sn we use some other effective action sffi in the same 
hierarchy we directly determine the associated estimator at this level. For 
example, the p = 2 level effective action equals [1,2] 



c(2) 



N-l 

E 

n=0 



5 



N 



+ eNV(q n ) + ^-V'( q - n ) + ^V"(q n ) 



2e N 



12 



24 



(13) 



Using the above procedure we easily determine the associated estimator to be 



E 



p=2 



1 

TV 



JV-l 

E 

n=0 



v + ^v' + —V" + -±V" + 

2 6 12 24 



'In , -/ < \ i ■// . <y „ i -II ( lr< f -N_ym _|_ Qn^n ym 



48 



, (14) 



where V is shorthand notation for V(q n ). For higher p levels one proceeds 
in precisely the same way. The effective actions become more complex as we 
move up the hierarchy, so that the corresponding p level estimators are most 
easily calculated using using a standard package for algebraic calculations such 
as MATHEMATICA. Explicit expressions for the energy estimators for p < 6 
are given in the Appendix. The results for p < 9 as well as the Monte Carlo 
codes used in this paper can be found on our web site [24]. 

The outlined procedure of matching a discretized action to its corresponding 
estimator has the property that amplitudes (and partition functions) have the 
same continuum limit as energy expectation values. As a result, by consistently 
using p level effective actions and associated estimators we are guaranteed to 
have 

U%\p) = U(P) + 0((P/Ny) . (15) 



In the following section we present explicit Monte Carlo simulations that verify 
this behavior. 



3 Numerical results 



The numerical simulations presented in this section were done using Grid- 
adapted Monte Carlo code and were run on EGEE-II and SEE-GRID-2 in- 
frastructure [25,26]. As already indicated, we carried out a series of Monte 
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Carlo simulations for several models, including anharmonic oscillator, mod- 
ified Poschl-Teller potential, and Morse potential. We determined explicitly 
that the speedup in convergence of energy expectation values of a general 
model precisely conforms to the analytical result given in Eq. (15). The fig- 
ures presented in this section give the results of simulations for the case of 
an anharmonic oscillator with quartic coupling using different levels p. The 
potential is thus 

V(q) = \q 2 + f g 4 ■ (16) 



All presented results correspond to g = 24 and [3=1, however, the same 
kind of increase in convergence is seen for all values of coupling and inverse 
temperature (and, in fact, for all models numerically analyzed). 
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Fig. 1. Discretized energy expectation values Uj^ as functions of N for p = 1,5,9 
for the case of an anharmonic oscillator with quartic coupling g = 24, time of prop- 
agation [3=1, and Nmc = 10 9 . The lines represent appropriate 1/N polynomial 
fits to the data. The level p curves have 1/N P leading behavior. 

(p) 

Fig. 1 shows how discretized energy expectation values U}^ converge to the 
continuum for p = 1 (naively discretized action and standard "virial" estima- 
tor), as well as for the effective actions and estimators at levels p = 5 and 
p = 9. From the figure we see that effective actions and estimators outper- 
form the naively discretized action and standard estimator. The lines represent 
curve fits of the data to polynomials in 1/N. We see that all p levels have the 
same continuum limit (within the error). The leading behavior of the p level 
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curve fit is 1/N P . The inset plot gives a zoomed in view clearly showing that 
for p = 9 one obtains excellent convergence (energy expectation value to five 
significant figures) for even extremely coarse discretizations such as N = 3. 

An even more clear-cut demonstration of the fact that numerical simulations 
conform to Eq. (15) can be seen in Fig 2. The plot display deviations from the 
continuum limit as a function of N. The dashed lines are polynomial data fits 
in 1/N starting from 1/N P , the leading terms 1/N P are given as solid lines. 
The deviations from the continuum limit — U\ become extremely small 
for larger values of N and p and are masked by statistical Monte Carlo errors. 
For this reason, one would need to use an even larger Nmc m order to show 
the p = 9 level deviation curve. Said another way, the p = 9 result for N = A 
is already well within the statistical error for Nmc — 10 9 . 



t 1 1 1 1 1 1 1 1 r 




2 4 8 16 32 64 128 256 512 1024 



Fig. 2. Log- log plot of the deviations from the continuum limit \U]E — U\ as functions 
of N for p = 1, 2 and 5 for an anharmonic oscillator with quartic coupling g = 24, 
time of propagation = 1, and Nmc = 10 9 . Dashed lines correspond to appropriate 
1/N polynomial fits to the data. Solid lines give the leading 1/N behavior showing 
that the level p curve has 1/N P leading behavior. 

The same behavior was seen for other values of parameters as well as for 
the two other models tested: particle in a modified Poschl- Teller potential 
Vmpriq) = — |[x 2 A(A — 1)/ cosh 2 (x<20] for x — 0-5 and A = 15.5, and particle 
in a Morse potential V M {q) = C{e- 2Kq - 2e- Kq ) for C = 10 and k = 2. 

As mentioned earlier, the effective actions and energy estimators for a give 
N get more complex as we go to higher values of p. Consequently, there is 
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a relative increase in computation time as we increase p. Fig. 3 shows this 
dependence as a function of p. We notice that for p > 5 the dependence is 
essentially exponential. Because of this increase in complexity, computations 
using Spf are about 37 times slower than that with S^. On the other hand, 
increase of p drastically improves convergence to the continuum limit, making 
it possible to obtain the same precision using much smaller values of N. Even 
for p = 9 the gain in precision outweighs the increase in computation time 
arising from complexity by many orders of magnitude. 



l i i i i i i i r 




123456789 



Fig. 3. Relative increase in computation time that comes from the increased com- 
plexity of expression for higher p-level effective actions and energy estimators for 
fixed N. For p > 5 this "complexity penalty" is well approximated by an exponen- 
tial - a direct consequence of the fact that p level effective actions follow from the 
starting action via a p-fold recursive process. 



To conclude, we have extended the method for systematically speeding up path 
integral calculation introduced in [1,2] to calculation of energy expectation 
values. We have shown that a consistent choice of effective action and estimator 
leads to the same form of speedup for expectation values that was previously 
seen for transition amplitudes, i.e. at p level the discretized expectation values 
differ from the continuum as 1/N P . 
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A Energy estimators to p = 6 level 



In the following we use the short hand notation q = q n , V = V(q n ), 5 2 = 5] 
and € — 6jy = P/N. 
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